# imports
import numpy as np
from aicspylibczi import CziFile
import matplotlib.pyplot as plt
import glob
import scipy.optimize
Background and Introduction¶
For our project we are using image analysis techniques from BIOL300 to obtain values from images of FRAP (fluorescence recovery after photobleaching) techniques over time. We can use these values along with techniques from BIOL301 to model the recovery curve of a specific FRAP dataset that we will choose. A part of our project will be modeling diffusion.
For this project, we chose to use Dr. Nikki Farnsworth's paper Modulation of Gap Junction Coupling Within the Islet of Langerhans During the Development of Type 1 Diabetes. This paper examines the role of gap junctions, specifically Connexin 36 (Cx36), in the role of Type 1 Diabetes (T1D). Therefore, the primary goal was to determine if changes in islet Cx36 gap junctions coupling and $Ca^{2+}$ signaling dynamics occur early in the progression of T1D, prior to changes in glucose tolerance. FRAP was used in this research paper so that one would be able to measure the rate at which fluorescently labeled molecules bound to Cx36 diffused between cells. As well as this, it was also used to help compare the FRAP recovery rates in islets from healthy and diabetic mice. The FRAP results help visualize how altering gap junction coupling can affect insulin secretion and $Ca^{2+}$ signaling.
A differential equation that models the behavior of F:
$$\frac{dC}{dt} = k_{off}(C_{tot} - C) + k_{on}C$$In this model, where C represents the concentration of fluorescent molecules in the bleached region at time t, $C_{tot}$ is the total concentration of fluorescent molecules, $k_{on}$ is the association rate constant, and $k_{off}$ is the dissociation rate constant. This equation will allow us to explore the exchange between bleached and unbleached molecules within our area of interest. $C_{tot}$ and values of C over time are the values that we can find using image analysis. Using the DE, we can then calculate the diffusion coefficients and create a recovery curve. Dr. Farnsworth's FRAP data that we are using can be seen in the Extracting Images subsection.
The intellectual merit of this project involves applying models with differential equations to the experimental technique of FRAP to graph and analyze FRAP-related data. FRAP can be a highly useful tool for examining movement in cells, including the diffusion of molecules, lateral and transverse movement of phospholipids in the cell membrane, and other processes. By modeling this graphically and mathematically, we will gain a better understanding of the data to be analyzed, potentially resulting in a better understanding in the biological process represented by the data.
The broader impact of this project is that FRAP analysis has the potential to make significant contributions to various fields including biological research, drug discovery, and environmental science. This project has the ability to provide insights into protein-protein interactions, cellular signaling, and enzyme activity. As well as this, it can help increase our understanding of processes such as diffusion, transportation, and more with cells. Next, FRAP analysis can help with identifying proteins involved in diseases by studying their mobility in cells. We can also use this information to see the impact of certain drugs on protein-protein interactions and cellular processes. Finally, FRAP analysis can be used to increase understanding on the impact of certain pollutants on cellular processes and organims health.
Defining Functions¶
extract_czi()¶
The purpose of the extract_czi() function is to process CZI image files and extract a sequence of image frames from them, each image frame representing a specific time point during the FRAP experiment. Then, these image frames will be organized into a list of image data by time.
def extract_czi(path):
'''Given a path to a czi file, reads data from the file into a list of images by time.'''
# read the czi file
czi = CziFile(path)
dims = czi.get_dims_shape()[0] # dictionary of dimensions (X, Y, and T are the ones we care about)
# get time dimension
time = dims['T']
# get the array of image data at each time point
images = []
for t in range(time[1]):
im = czi.read_image(T=t)
images.append(im[0])
# remove extra dimensions of each image (squeeze)
for i in range(len(images)):
images[i] = np.squeeze(images[i])
return images
# see Image Analysis subsection for testing
graph_subplots()¶
The purpose of the graph_subplots() is to visualize a list of images as a grid of subplots. This function will enable us to see the time series data.
def graph_subplots(data, num_rows, num_columns, title=True):
'''Takes image data and graphs it as subplots.'''
# set axes
fig, axs = plt.subplots(num_rows, num_columns, figsize=(15,15))
plt.tight_layout()
# graph each image
for i, ax in enumerate(axs.flat):
if i < len(data):
ax.imshow(data[i])
# label with time step if desired
if title:
ax.set_title(f"Time Step t={i+1}")
# see simulate diffusion test case for testing
slice_image()¶
The purpose of the slice_image() function is to help extract a rectangular portion from a given image.
def slice_image(image, x, y, width, height):
'''Slices and returns a rectangular portion of the given image.'''
# initialize list to store sliced image data
sliced = []
# add each pixel in the region to slice
for i in range(y, y+height):
new_row = []
for j in range(x, x+width):
new_row.append(image[i][j])
sliced.append(new_row)
return np.array(sliced)
# see Image Analysis subsection for testing
average_fluorescence()¶
Purpose
def average_fluorescence(images):
'''Given a list of images, returns a list of the average fluorescence for each.'''
# initialize list to store average fluorescences
avg_fluorescences = []
# iterate through each image
for im in images:
im_fluorescence = []
# store the average fluorescence of each row
for row in im:
im_fluorescence.append(np.mean(row))
# find the total average fluorescence of the image and add it to the list of average fluorescences
avg_fluorescences.append(np.mean(im_fluorescence))
return avg_fluorescences
# see Comparison of the Average Fluorescence Over Time Between the Unbleached vs Bleached Regions for testing
find_nearest_pixel()¶
The purpose of the find_nearest_pixel() is to help find the nearest fluorescence pixel in a specific direction from the current pixel.
def find_nearest_pixel(image, location, direction, fluorescable_map):
'''Given an image and the location of a pixel in that image, returns tuple (x,y) of the
location of the nearest fluorescable pixel within an 11x11 block in a given direction,
or -1 if a pixel cannot be found.'''
# extract x and y coordinates from location
x_loc = location[0]
y_loc = location[1]
# build the 11x11 block to check and a corresponding fluorescence map
# initialize lists to store the block and fluorescence map
block = []
block_map = []
# find indices of the edges of the block
left_edge = x_loc - 5
if left_edge < 0: # in case block hits the bounds of the image
left_edge = 0
right_edge = x_loc + 6
if right_edge > len(image[0]):
right_edge = len(image[0])
up_edge = y_loc - 5
if up_edge < 0:
up_edge = 0
down_edge = y_loc + 6
if down_edge > len(image):
down_edge = len(image)
# create the block by adding the pixels from the original image
y_count = 0
for y in range(up_edge, down_edge):
block_row = []
map_row = []
x_count = 0
for x in range(left_edge, right_edge):
# store the new location of the pixel we're searching from
if x == x_loc and y == y_loc:
block_x = x_count
block_y = y_count
# add the corresponding pixel from the image to the block and the map
block_row.append(image[y][x])
map_row.append(fluorescable_map[y][x])
x_count += 1
block.append(block_row)
block_map.append(map_row)
y_count += 1
# initialize variables to store the minimum directional distance, minimum total distance, and nearest pixel location
min_dist = np.inf
min_total = np.inf
nearest_p = -1
# find nearest pixel within block
for y, row in enumerate(block):
for x, p in enumerate(row):
# skip any pixels that aren't fluorescable based on the fluorescence map
if not block_map[y][x]:
continue
# calculate directional and total distances
left_dist = block_x - x
right_dist = -left_dist
up_dist = block_y - y
down_dist = -up_dist
total_dist = (left_dist**2 + up_dist**2)**0.5
# depending on selected direction, stores the location of the pixel if it's closer than the existing closest
if direction == "left":
if left_dist <= min_dist and left_dist > 0 and total_dist <= min_total:
nearest_p = x, y
min_dist = left_dist
min_total = total_dist
# catch left border
if location[0] == 0:
return -1
elif direction == "right":
if right_dist <= min_dist and right_dist > 0 and total_dist <= min_total:
nearest_p = x, y
min_dist = right_dist
min_total = total_dist
# catch right border
if location[0] == len(image[0]):
return -1
elif direction == "up":
if up_dist <= min_dist and up_dist > 0 and total_dist <= min_total:
nearest_p = x, y
min_dist = up_dist
min_total = total_dist
# catch top border
if location[1] == 0:
return -1
elif direction == "down":
if down_dist <= min_dist and down_dist > 0 and total_dist <= min_total:
nearest_p = x, y
min_dist = down_dist
min_total = total_dist
# catch bottom border
if location[1] == len(image):
return -1
# convert coordinates of nearest pixel to be relative to the image rather than the block
if nearest_p != -1:
x_dist = nearest_p[0] - block_x
y_dist = nearest_p[1] - block_y
nearest_x = x_loc + x_dist
nearest_y = y_loc + y_dist
nearest_p = nearest_x, nearest_y
return nearest_p
# test the block construction of find nearest pixel function
def test_check_block(x_loc, y_loc, image):
# build the 11x11 block to check and a corresponding fluorescence map
# initialize lists to store the block and fluorescence map
block = []
block_map = []
# find indices of the edges of the block
left_edge = x_loc - 5
if left_edge < 0: # in case block hits the bounds of the image
left_edge = 0
right_edge = x_loc + 6
if right_edge > len(image[0]):
right_edge = len(image[0])
up_edge = y_loc - 5
if up_edge < 0:
up_edge = 0
down_edge = y_loc + 6
if down_edge > len(image):
down_edge = len(image)
# create the block by adding the pixels from the original image
y_count = 0
for y in range(up_edge, down_edge):
block_row = []
# map_row = [] ignoring map for test case
x_count = 0
for x in range(left_edge, right_edge):
# store the new location of the pixel we're searching from
if x == x_loc and y == y_loc:
block_x = x_count
block_y = y_count
# add the corresponding pixel from the image to the block and the map
block_row.append(image[y][x])
# map_row.append(fluorescable_map[y][x])
x_count += 1
block.append(block_row)
# block_map.append(map_row)
y_count += 1
return block
image = [
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
]
x_loc = 6
y_loc = 7
test_check_block(x_loc, y_loc, image)
[[1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1], [1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0], [1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]]
# test find nearest pixel function
test_image = [
[0, 0, 0, 1, 1],
[0, 0, 0, 0, 0],
[0, 0, 1, 0, 0],
[0, 0, 0, 0, 0],
[1, 0, 0, 1, 0]
]
test_map = [
[0, 0, 0, 1, 1],
[0, 0, 0, 0, 0],
[0, 1, 1, 0, 0],
[0, 0, 0, 0, 0],
[1, 0, 0, 1, 0]
]
test_loc = (2, 0) # 0 in middle of first row
print(f"Left: {find_nearest_pixel(test_image, test_loc, 'left', test_map)}") # expected (1, 2)
print(f"Right: {find_nearest_pixel(test_image, test_loc, 'right', test_map)}") # expected (3, 0)
print(f"Up: {find_nearest_pixel(test_image, test_loc, 'up', test_map)}") # expected -1
print(f"Down: {find_nearest_pixel(test_image, test_loc, 'down', test_map)}") # expected (2, 2)
Left: (1, 2) Right: (3, 0) Up: -1 Down: (2, 2)
simulate_diffusion()¶
The purpose of simulate_diffusion() is to help simulate the process over time given an initial image and fluorescence map. This will help us utilize a probabilistic approach to model the movement of the fluorescent molecules between pixels.
def simulate_diffusion(image, fluorescable_map, delta_t, num_steps, k):
'''Given a starting image and map of fluorescable pixels, simulates diffusion using probability and returns
array of image data to display the results of the simulation.'''
# get height and width of image
height = len(image)
width = len(image[0])
# build and initialize probability array
p = np.zeros((num_steps, height, width))
# initial condition is the fluorescence of our starting image--more fluorescent means it's more likely to diffuse
for y in range(height):
for x in range(width):
p[0][y][x] = image[y][x]
# model diffusion over time
for t in range(num_steps-1):
# loop over each pixel
for y in range(height):
for x in range(width):
# don't diffuse if the pixel can't fluoresce
if fluorescable_map[y][x] == 0:
continue
pix_loc = (x, y)
# find the nearest pixel in each direction--these are the "boxes" fluorescence can diffuse to
left = find_nearest_pixel(image, pix_loc, "left", fluorescable_map)
right = find_nearest_pixel(image, pix_loc, "right", fluorescable_map)
up = find_nearest_pixel(image, pix_loc, "up", fluorescable_map)
down = find_nearest_pixel(image, pix_loc, "down", fluorescable_map)
# store each of the valid pixels
boxes = [left, right, up, down]
diffusable_boxes = [box for box in boxes if box != -1]
# initialize new probability
new_prob = 0
# add to the new probability for fluorescence diffusing into our current pixel
for box in diffusable_boxes:
new_prob += p[t, box[1], box[0]]
# remove from the new probability for fluorescence diffusing out of our current pixel
new_prob -= len(diffusable_boxes) * p[t, y, x]
# update probability matrix
p[t+1, y, x] = p[t, y, x] + k*delta_t*new_prob
return p
# test diffusion simulation
t_steps = 20
dt = 0.1
k_guess = 1
test_image = [
[0, 1, 0, 1, 1, 0, 1, 1, 1, 0, 0],
[1, 1, 1, 1, 1, 1, 1, 0, 1, 0, 0],
[1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 1],
[0, 1, 1, 0, 1, 1, 1, 0, 0, 1, 1],
[1, 0, 0, 1, 1, 1, 1, 1, 0, 1, 1],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0],
[0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0]
]
test_map = [
[0, 1, 0, 1, 1, 0, 1, 1, 1, 0, 0],
[1, 1, 1, 1, 1, 1, 1, 0, 1, 0, 0],
[1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 1],
[0, 1, 1, 0, 1, 1, 1, 0, 0, 1, 1],
[1, 0, 0, 1, 1, 1, 1, 1, 0, 1, 1],
[0, 1, 1, 1, 1, 0, 1, 1, 0, 1, 0],
[1, 1, 1, 0, 1, 1, 0, 0, 1, 1, 1],
[1, 0, 1, 0, 1, 1, 0, 1, 0, 0, 1],
[0, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0]
]
p = simulate_diffusion(test_image, test_map, dt, t_steps, k_guess)
graph_subplots(p, 5, 4)
CZI class¶
The purpose of the CZI class is to help with future applications of this project. By using the set_bleached_dims() and set_unbleached_dims() functions allows one to define separate bounding boxes for bleached and unbleached regions within the chosen (FRAP czi) image. Then, the slice_image() function iterates through each image and slices it based on the defined region. This stores the unbleached, bleached, and combined data in separate lists.
class CZI:
def __init__(self, path):
'''Constructs CZI object using file path'''
# store critical information
self.path = path
self.images = extract_czi(self.path)
self.time_steps = len(self.images)
# initialize dimensions for bleached and unbleached slices
self.x_b = None
self.y_b = None
self.width_b = None
self.height_b = None
self.x_u = None
self.y_u = None
self.width_u = None
self.height_u = None
self.bleach_ind = None
# initialize lists for sliced regions
self.sliced_data = []
self.unbleached_data = []
self.bleached_data = []
def set_bleach_ind(self, bleach_ind):
'''Sets and stores the time index at which photobleaching occurred.'''
self.bleach_ind = bleach_ind
def set_bleached_dims(self, x, y, width, height):
'''Sets and stores the dimensions and location of the bleached region of the images.'''
self.x_b = x
self.y_b = y
self.width_b = width
self.height_b = height
def set_unbleached_dims(self, x, y, width, height):
'''Sets and stores the dimensions and location of the unbleached region of the images.'''
self.x_u = x
self.y_u = y
self.width_u = width
self.height_u = height
def slice_images(self):
'''Uses bleached and unbleached dimensions to slice the image data and store it in the object.'''
# loop over each image and slice it
for im in self.images:
unbleached_slice = slice_image(im, x_u, y_u, width_u, height_u)
bleached_slice = slice_image(im, x_b, y_b, width_b, height_b)
# store sliced data
self.unbleached_data.append(unbleached_slice)
self.bleached_data.append(bleached_slice)
self.sliced_data.append(np.concatenate((unbleached_slice, bleached_slice)))
def show_images(self):
'''Displays every image in the czi file.'''
to_graph = []
for image in self.images:
to_graph.append(image.data)
graph_subplots(to_graph, 5, 5)
def __str__(self):
'''Establishes string to display when object is printed.'''
return self.path
Image Analysis¶
To conduct image analysis we began by inputting Dr. Farnsworth's FRAP czi data from her study Modulation of Gap Junction Coupling Within the Islet of Langerhans During the Development of Type 1 Diabetes. Once the data was input, we then extracted the image data from every czi in the folder. To make sure our code was working properly we tested the first image in the folder, FRAP_data/02102017\RIP_0.01_FRAP1.czi, where RIP stands for Rat Insulin Promoter. It is important to note that the following code can easily be applied to any of the images in the folder, the slices and indices would have to be adjusted.
# get files
path_list = glob.glob("FRAP_data/02102017/*.czi")
path_list
['FRAP_data/02102017\\RIP_0.01_FRAP1.czi', 'FRAP_data/02102017\\RIP_0.01_FRAP2.czi', 'FRAP_data/02102017\\RIP_0.01_FRAP3.czi', 'FRAP_data/02102017\\RIP_0.1_FRAP1.czi', 'FRAP_data/02102017\\RIP_0.1_FRAP2.czi', 'FRAP_data/02102017\\RIP_0.1_FRAP3.czi', 'FRAP_data/02102017\\RIP_0_FRAP1.czi', 'FRAP_data/02102017\\RIP_0_FRAP2.czi', 'FRAP_data/02102017\\RIP_0_FRAP3.czi', 'FRAP_data/02102017\\WT_0.01_FRAP1.czi', 'FRAP_data/02102017\\WT_0.1_FRAP1.czi', 'FRAP_data/02102017\\WT_0.1_FRAP2.czi', 'FRAP_data/02102017\\WT_0_FRAP1.czi', 'FRAP_data/02102017\\WT_0_FRAP2.czi', 'FRAP_data/02102017\\WT_0_FRAP3.czi']
# store each file as a czi object
czi_files = []
for path in path_list:
czi_files.append(CZI(path))
# select czi file to work with and display
czi = czi_files[0]
czi.show_images()
The above images are our dataset which captures the fluorescence intensity over time. The images will then be sliced into two regions: bleached (which will initially be dark after photobleaching) and unbleached (which will retain the fluorescence).
# define coordinates for unbleached region to slice
x_u = 98
y_u = 82
width_u = 307
height_u = 180
czi.set_unbleached_dims(x_u, y_u, width_u, height_u)
# define coordinates for bleached region to slice
x_b = 98
y_b = 262
width_b = 307
height_b = 215
czi.set_bleached_dims(x_b, y_b, width_b, height_b)
# slice the images and store
czi.slice_images()
unbleached_data = czi.unbleached_data
bleached_data = czi.bleached_data
sliced_data = czi.sliced_data
# display sliced images
graph_subplots(sliced_data, 5, 5)
The above sliced images are a more helpful way to show fluorescence recovery in the Rat Insulin Promoter. As we can see in the above images, they have been sliced into two regions: bleached and unbleached. The top half of the images have retained fluorescence meaning that this is the unbleached region. Then, the bottom half which appears darker as the time step increases is darker, meaning that it is the bleached region.
To zoom in on the meaning of some of these images we can examine some of the time steps: t=3, t=4, and t=22. At the t=3 time step, the image captured shows a pre-bleached fluorescence. Then, at t=4 the image captures what the RIP looks like right after photobleaching. We can see that the bleached region decreases significantly. Finally, at the t=22 time step we can see that fluorescence is beginning to recover in the bleached region, this is caused by diffusion.
# show unbleached, bleached, and total sliced regions immediately after photobleaching has occurred
czi.set_bleach_ind(3)
bleach_ind = czi.bleach_ind
slice_zooms = [unbleached_data[bleach_ind], bleached_data[bleach_ind], sliced_data[bleach_ind]]
graph_subplots(slice_zooms, 1, 3, title=False)
The above images are zooming in on the slices to show that slicing was successful. We can see that our top half is unbleached because fluorescence was retained and that the bottom half is bleached because it appears darker. The increase in fluorescence intensity in the bleached region is consistent with the expected recovery due to molecular diffusion. As well as this, the decrease in fluorescence intensity in the unbleached region can correspond to molecules diffusing out to restore the equilibrium. Finally, these images are showcasing that molecular diffusion is occuring as expected.
Comparison of the Average Fluorescence Over Time Between the Unbleached vs Bleached Regions¶
# calculate average fluorescences of unbleached and bleached regions
unbleached_fluorescences = average_fluorescence(unbleached_data)
bleached_fluorescences = average_fluorescence(bleached_data)
# get time values
t = np.arange(czi.time_steps)
# graph average fluorescence over time for each region
plt.plot(t, bleached_fluorescences, '.', label="Bleached Region")
plt.plot(t, unbleached_fluorescences, '.', label="Unbleached Region")
plt.xlabel("Time")
plt.ylabel("Average Fluorescence")
plt.legend(loc="right")
plt.show()
In the above graph, the average fluorescence was calculated for both the unbleached and bleached regions across all of the time steps. From the comparison of the average fluorescence over time we can see that in the unbleached region decreases slightly over time as the molecules diffuse out. In the bleached region, we can see that it starts to dwindle, then rapidly decreases, but as time continues it is recovered. Therefore, this graph aligns with the expected behavior of FRAP experiments, which helps to confirm the diffusion of fluorescent molecules from the unbleached region to the bleached region over time.
Simulated Diffusion¶
One of the ways we modeled FRAP was with a probability-based simulation. The concept is derived from the box-jumping master equation model, where we create a matrix of the probability a particle is in a given “box” and update that matrix over a set number of time steps. We expanded this model from a linear one-dimensional model to a two-dimensional model so we could use it to model the diffusion in this two-dimensional cell.
We took the first image after photobleaching occurred and thresholded it to have binary fluorescence values (0 or 1), then used that as the start image for our simulation. The simulation then, for each pixel in the image, finds the nearest “fluorescable” pixel above, below, to the left, and to the right of that pixel. These pixels are determined to be fluorescable based on a fluorescence map, which is another thresholded image from the start of the experiment. The fluorescence of each of these pixels is directly related to how likely it is to diffuse, so the simulation uses probability to adjust the fluorescence of each. It does this for every time step, then we can examine the grid of images it generates to see the simulation of diffusion over time.
# create a fluorescability map by thresholding the first image in the file
map_base = sliced_data[0]
thresh = 4000
fluorescable_map = map_base > thresh
plt.imshow(fluorescable_map);
# create an image to start the diffusion simulation by thresholding the image immediately after photobleaching has occurred
start = sliced_data[bleach_ind]
thresh = 4000 # same threshold as map for consistency
sim_start = start > thresh
plt.imshow(sim_start);
# set parameters
t_steps = czi.time_steps - bleach_ind # run simulation for as many time steps as the czi has after photobleaching
dt = 0.1
k_guess = 1
# run diffusion simulation and display results
p = simulate_diffusion(sim_start, fluorescable_map, dt, t_steps, k_guess)
graph_subplots(p, 5, 4)
# compare beginning and end of diffusion simulation side by side
to_compare = (p[0], p[-1])
graph_subplots(to_compare, 1, 2, title=False)
In the images above, which are the beginning and end of the simulation respectively, we can see that fluorescence has spread out through the simulated cell to start filling gaps, as would be expected in the diffusion process. These can then be compared to the beginning and end of diffusion in the original image data, which is shown below:
# show beginning and end of original data side by side
to_compare = [sliced_data[bleach_ind], sliced_data[-1]]
graph_subplots(to_compare, 1, 2, title=False)
Differential Equations¶
The goal is to make a plot of the integrated differential equation $\frac{dC}{dt} = k_{off}(C_{tot} - C) + k_{on}C$, which describes FRAP behavior, use this as our predicted curve, and then compare it to our actual curve. $C$ will be our fluorescence values over time, $C_{tot}$ will be the total concentration of fluorescent molecules, $k_{off}$ will be the dissociation rate constant, and $k_{on}$ will be the association rate constant.
First, we have defined a few functions to help us plot both our predicted curve and our actual curve. The average_fluorescence function will give us a list of average fluorescent values from the given image files. This data can plotted as our actual curve over time. The FRAP_rhs function defines the differential equation that we will be integrating to find our predicted FRAP curve.
SORRY -- I moved average fluorescence up to the function definitions so this may need to change a little bit
def FRAP_rhs(C, t, C_tot, k_off, k_on):
'''Describes right-hand side of FRAP differential equation.'''
# compute dC/dt
dC_dt = k_off*(C_tot - C) - k_on*C
# Return the result as a NumPy array
return np.array(dC_dt)
After defining our functions, we can plot our predicted and actual data over time on the same graph. This will help us visualize the behavior of our actual FRAP data versus what we would expect it to do. We first defined the initial values and parameters of the differential equation we are going to integrate. For now, some of these will be guesses that we can later find optimized values for. Then, using the scipy.integrate.odeint function, we integrated our DE to create our predicted curve (C_guess). The predicted curve and the actual data are plotted against time on the same plot below. We can see visually that the predicted curve fits our data fairly well.
# integrate diff eq and graph with guesses
# get data (average fluorescences after bleaching) and time points
data = np.array(average_fluorescence(czi.images[bleach_ind:]))
t = np.arange(len(data))
# initial condition guess
C_0_guess = data[0]
# parameter guesses
C_tot_guess = data[0]
k_off_guess = 0.15
k_on_guess = -0.008
# package parameters into a tuple
guess_args = (C_tot_guess, k_off_guess, k_on_guess)
# integrate diff eq
guess_fit = scipy.integrate.odeint(FRAP_rhs, C_0_guess, t, args=guess_args)
# compare actual data to predicted curve found by numerical integration
plt.plot(t, data, '.', label="Data")
plt.plot(t, guess_fit, label="Guess Fit")
plt.xlabel("Time")
plt.ylabel("Fluorescence")
plt.legend()
plt.show()
After doing this with "best guess" values for some of our parameters, we then found optimized values. Below we defined a function to return the chi-squared value for our provided parameters and data. We then used this function to perform scipy.optimize.minimize in order to find the optimal values we should use based on our guesses.
def chi_squared_FRAP(params, data, t):
'''Returns the chi squared value for the parameters and data provided.'''
# unpack the parameters (guesses)
C_tot = params[0]
k_off = params[1]
k_on = params[2]
C_0 = params[3]
arguments = (C_tot, k_off, k_on)
# compute the fit by integrating diff eq
fit = scipy.integrate.odeint(FRAP_rhs, C_0, t, args=arguments)
# compute the chi-squared
vals = (data - fit)**2 / fit
return np.sum(vals)
# optimize values by minimizing chi squared
optimal_vals = scipy.optimize.minimize(chi_squared_FRAP, [C_tot_guess, k_off_guess, k_on_guess, C_0_guess], args=(data, t),
bounds=[ (2000, 2500), (-1, 1), (-1, 1), (2000, 2500) ])
# extract optimal values
best_C_tot = optimal_vals.x[0]
best_k_off = optimal_vals.x[1]
best_k_on = optimal_vals.x[2]
best_C_0 = optimal_vals.x[3]
print(f'Best C tot: {best_C_tot}, Best k off: {best_k_off}, Best k on: {best_k_on}, Best C_0: {best_C_0}')
Best C tot: 2006.3229817087017, Best k off: 0.9999902430202654, Best k on: -0.0712297307784808, Best C_0: 2160.1942355373376
After finding the optimal values, we plugged them back into our DE and once again integrated it to find the curve of best fit. We plotted this curve over time alongside our actual data and our initial guess fit curve. As seen below, this plot shows how the best fit compares to the guess fit and matches the actual data better. This method demonstrates how a FRAP curve can be predicted using an optimization and integration of this particular FRAP differential equation model.
# calculate best fit by integrating with optimized values
best_args = (best_C_tot, best_k_off, best_k_on)
best_fit = scipy.integrate.odeint(FRAP_rhs, best_C_0, t, args=best_args)
# calculate polynomial fit
coeffs = np.polyfit(t, data, 2)
poly_func = np.poly1d(coeffs)
poly_fit = poly_func(t)
# plot data and fit curves
plt.plot(t, data, '.', label="Data")
plt.plot(t, guess_fit, label="Guess Fit")
plt.plot(t, best_fit, label="Scipy Fit")
plt.plot(t, poly_fit, label="Poly Fit")
# label and show graph
plt.xlabel('Time')
plt.ylabel('Fluorescence')
plt.legend()
plt.show()
Results and Discussion¶
Image analysis, diffusion simulation, graphing. Will do this when all code is up and running correctly, we can talk about why each thing helped us analyze the FRAP data or why it is a useful model.
Next Steps¶
If we were to explore further methods of analyzing FRAP data, we could start by comparing the FRAP curves of different parts of the cell to each other as well as curves from the same part of the cell under different conditions. For instance, the dataset we used to for this project contained multiple files of FRAP data that were measured under different conditions. The purpose was the track the diffusion of Cx36 between the cells of mice that were genetically edited to display different diabetic conditions. We could use the methods of modeling that we developed to predict the behavior of Cx36 diffusion under each set of conditions. When comparing this to the actual data, the results should match closely to support whatever effect is observed.
Need to continue this